Uniform Approximation of Linear Systems* 
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A method for reducing the complexity of the class of linear, time-varying, 
dynamic control systems is developed where the approach taken is that of 
uniform approximation (that is, modeling for a region of initial conditions). 
The objective of the modeling procedure is to choose a linear time-invariant 
system of given dimension, that minimizes a "worst-case" type of error 
criterion. Some results from the theory of widths of sets in Banach space 
are used to obtain bounds on the optimal approximation error as a function 
of the dimension of the approximating system. The use of these bounds in 
choosing the order of the approximation is discussed. An example illustrates 
the use of the derived results. 

I. INTRODUCTION 

In the analysis and design of control systems it is often useful to 
have low order constant coefficient models for the system. The prob- 
lem of modeling linear systems by lower order linear systems has 
received considerable attention, but these analyses have usually been 
restricted to the modeling of constant coefficient systems. 

References 1 through 5 contain various approaches to the system 
approximation problem; however, these analyses are generally re- 
stricted to the modeling of constant coefficient systems or systems 
which are forced with a given input or initial condition. 

The control system analyst often finds himself dealing with non- 
stationary systems, but little work has been done in the area of 
optimally modeling this class of systems. The emphasis here is on 
modeling the class of linear, homogeneous time-varying systems with 
constant coefficient models. Reference 6 considers approximation of 
forced systems. Rather than design the model requiring solutions of 
the actual and approximate systems be "close" for a prescribed initial 
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condition, the approach taken here is that of uniform approximation. 
Initial conditions are assumed to lie in some set in Euclidean space 
and a "worst-case" type of error criterion is denned. This eliminates 
tuning the model to specific conditions which may not be met when 
using the model. The material presented here thus generalizes previous 
work in that it extends the class of systems considered to time-varying 
systems and generalizes the error criterion to handle the more realistic 
problem of modeling for regions of initial conditions. 

The problem is of importance, for example, in trajectory analysis 
where the linear time-varying system is obtained by linearizing a set 
of nonlinear equations about a nominal trajectory. In this case the 
time-varying nature of the system is described by partial derivatives 
evaluated along the nominal trajectory. Solutions to the resulting 
equations require simulation for each set of initial conditions. Using 
a constant coefficient model eliminates the need for repeated simula- 
tion. 

The above example illustrates the use of a simplified model in 
analysis. The designer is interested in reducing the complexity of 
high-order nonstationary control system plants since this provides 
a means for designing simpler controllers based upon the model de- 
scription. The results presented here not only allow one to obtain 
stationary models but simultaneously offer the opportunity to obtain 
lower order models of the original system. 

II. PROBLEM DEFINITION AND FORMULATION 

The system we are considering is described by the linear, time- 
varying, homogeneous vector differential equation 

x(t) = A(t)x(t) (1) 

with the outputs given by 

2/(0 = C(t)x(t) (2) 

where 

x(f) is an n-vector 

A(t) is an n X n matrix whose elements are bounded and piecewise 
continuous on [t , t f ]. 

C(t) is an m X n matrix whose elements are bounded and piecewise 
continuous on [t , t,\. 
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It is desired to obtain a constant coefficient system of /cth order* 
(fc ^ m) 

£(t) = Az{t) (3) 

such that the first m components of the state vector $(t) closely approxi- 
mate the components of y(t) over the finite time interval [t , t r ). Writing 

yit) = C£(t) (4) 

with 

C = [I mXm i 0] 

the approximation problem can be viewed as choosing the elements of 
the fc X k matrix A such that y(t) approximates y(t) over [t , t f ]. 

Since, in general, it is not known at the time of modeling what initial 
conditions will exist in the system, it is desirable to have the approxi- 
mating system depend on a prescribed range of initial conditions rather 
than being tuned to any specific initial condition. The initial conditions 
are considered to he in a closed, bounded convex subset of Euclidean 
n-space. That is, 

x(t ) t R C E n , 

and the performance criterion is given by 

J k (l) = max min [ ' {y - y)'W(t)(y - y) dt (5)t 

where 

[t , tf] is bounded 

y(t) is the solution of (1) and (2) with x(t ) = x 
y(t) is the solution of (3) and (4) with Z(t ) = £„ 
W{t) is positive definite and bounded for all t z [t , t f ]. 

The above performance criterion corresponds to the worst case error 
in the approximation, corresponding to a given model, when the initial 
condition on the model, $(t ), is chosen optimally in terms of the initial 
conditions on the actual system. The modeling objective is to choose 
A to minimize J k (A) (that is, minimize the maximum approximation 
error) . 
The approximation problem will be cast into a Hilbert space setting 



* Notice that k is not restricted from above. It may be desirable to have k 
> n if the original system is time-varying. 
t In all that follows the prime denotes transpose. 
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which will permit the use of many of the general results to be presented 
in the next section. Vector spaces of solutions of the original system 
equations and any member of the class of approximate system equations 
are established. These spaces are then imbedded into an encompassing 
Hilbert space. It then is shown that the problem of finding an optimal 
approximation can be viewed as a problem of finding the "best" subspace 
(of a given form) of the Hilbert space to use in approximating solutions 
of the original system. Writing the output vector of the original system 
in terms of the transition matrix leads to 

y(t) = emit, 0*(0 (6) 

where the transition matrix $(t, t ) satisfies 

| *<*, = A(t)Ht, Q (7) 

with initial conditions 

*(*. , Q - /. (8) 

Now if the original system is completely observable 8,9 on the finite 
interval [t , t f ] the columns of the m X n matrix C(t)$(t, t ) are linearly 
independent as vector-valued time functions. That is, 

C(t)$(t, t )x(t ) = for all 1 1 [t , t f ] 

implies x(t a ) = 0. For an observable system, the initial state can be 
determined uniquely from knowledge of the output. Since x(t ) = =► 
y(t) = and, from observability, y{t) = =» x(t ) = the linear in- 
dependence of the columns of C(t)$(t, t ) follows. 

Let <y be the linear space spanned by the n columns of C(t)$(t, t ). 
The solutions of the original system he in % which is of dimension n for 
an observable system. Notice that the number of components (w) in 
the vector y and the dimension of the space ^ need not be the sama If 
the system is not completely observable on [t , t f ] the dimension of °y is 
less than n. 

The solutions to equations (3) and (4) can be written as 

y(t) = Ce^'-'-'xiQ (9) 

where 

e ~ *-> VI 

t-0 * ! 

and 
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A^c-..) = Xe 2lt -"\ (10) 

at 

It is thus seen that solutions p(t) lie in the vector space spanned by the 
k columns of the in X k matrix cV <,_ ' 0> . Denote this vector space as 
% . If X is such that the approximate system is observable then «y* is 
of dimension k and the k columns of cV (, ~' o) form a basis 

{g. ;i= 1, ••• ,k) 

for the fc-dimensional vector space <y A of approximating solutions. These 
basis elements can be written as 

*(0 = Ce I( '-' 0, K,- (11) 

g, = MO;* e [«.,*/]) (12) 

where K { is the ith column of the k X fc identity matrix. If the approxi- 
mation is not observable the dimension of <y k is less than k. In any case 
vector spaces <y k with basis elements of the form (11) characterize the 
approximating systems where A is a k X k real matrix. Denning 

©* = fy*;gi, •'• ,g*span«y*} (13) 

where p<(<) is given by equation (11) and 1 is any real constant k X k 
matrix casts the problem into finding an element of £>* minimizing J k . 

The problem of finding an optimal approximation has been cast into 
the problem of finding an extremal space ^% e D* of approximating 
solutions. A Hilbert space 3C containing <y and all members of 3D* will 
now be constructed. _ 

Recall that the elements of «y and <y fc are real, vector-valued, time 
functions having m components. Thus each element of the Hilbert 
space 3C to be constructed will have m components. The inner product in 
3C is defined by 

(y, , y 2 ) = [ " yi(t)W(t)y a (t) dt (14) 

where W(t) is a real symmetric mX m matrix which is positive definite 
for t e [t„ , t f ] and whose elements are bounded for tt[t of t f ]. Notice that 
this is the same matrix appearing in the performance criterion given 
by equation (5). The norm of an element in 3C is given by 

llyll = (y, y)*. ( 15 ) 

The Hilbert space 3C is defined as 

3C = {y; y has m components, || y || < «» } 
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where || y || is given by (15) and the inner product given by (14). 
Since 

t f - L < 00 

and the elements of A(t) and C(t) are bounded it follows that solutions 
of equations (1) and (2) are bounded thus yielding 

Since elements of ^ k are bounded over the finite interval [t , t f ] 

That ^k and <y are subspaces of 3C follows from the fact that any finite- 
dimensional linear set in a normed space is closed 10 . 

The set of functions to be approximated are solutions to the original 
system equations with the initial conditions x(t„) satisfying 

x{Q tRCEn 

where R is a closed, bounded convex subset of Euclidean n-space. 
Writing 

* = {y; y(t) = C(0*(«, Qx(t ), x(t.) * R\ (16) 

gives 

J k (l) = max min || y - y || 2 (17) 

where the modeling objective is to find 

dl = inf max min || y - y || 2 . (18) 

Before proceding to solve the formulated approximation problem, 
some results from the theory of widths in Banach space are outlined. 
Lower bounds on the optimal performance are found as a function of 
the dimension of the approximating system. 

III. WIDTHS OF SETS IN BANACH SPACE AND LOWER BOUNDS* 

Classically, approximation theory was concerned with the follow- 
ing problem. Given a function to approximate and a set of approxi- 
mating functions (sinusoids, exponentials, and polynomials, for ex- 
ample) find that linear combination of approximating functions which 



* Ref . 7 contains an excellent treatment of widths of sets in Banach space. 
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minimizes some distance function. Notice that here the approximat- 
ing functions are given as part of the problem statement. 

Ra,ther than approximate a single function, the problem under con- 
sideration is to approximate the class of functions tf given by (16). 
For a given class of functions JF it is desired to obtain a "best" set of 
approximating functions rather than to choose the set arbitrarily. A 
measure of comparison is introduced which enables one to evaluate the 
efficiency of different sets of approximating functions. The following 
definitions serve to illustrate these ideas. 

Let (B be a Banach space containing a set of functions 5 to be approxi- 
mated by elements of an n-dimensional subspace, X n , of CB. It is desired 
to find the "best" n-dimensional subspace, or equivalently the "best" 
set of approximating functions to use in approximating elements of 3\ 

For a given / t 5 and X n C ® 

inf ||/ -x || 

xtX» 

represents how well one can do in approximating a given / with elements 
of X n . Taking the supremum of the above quantity over all elements 
in "5 leads to the following definition. 

Definition 1: The deviation of ff from X n is given by 

E X n(&) = sup inf || / — x \\. 

ftS xtX n 

The deviation represents the worst case approximation error over the 
class ff when using elements of X n . Notice that the deviation serves as 
a performance measure of X„ . Taking the infimum of the deviation 
over all n-dimensional subspaces of (B leads to the following definition. 

Definition 2: The nth width of SF is given by 
d n ((F) = inf E Xn ($) 

X„c(B 

= inf sup inf || / — x ||. 

.Y„C(B ftS XtXn 

Some of the elementary results following from the above definitions are 

(i) The mono tonicity of the width: 

do(5F) ^ di(*) ^ d *&) ^ ••• 
and 
(ii) The nested property: If ff, C ^2 C • • • then 

4(ffO ^ d n (Z 2 ) S ■■ ■ 
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Notice that 

-/*(!) - El k (S). (19) 

In defining d\ the infimum of the square of the deviation was taken 
over the j'-dimensional (j ^ fc) subspaces in SD* whereas in denning d k 
the infimum was taken over all A;-dimensional subspaces of (B. Using the 
monotonicity property of the width, with 3C serving as the required 
Banach space, gives 

J t (l) £ df(ff) (20) 

for any k X k matrix A. 

Definition 3: U„ is a closed ball of radius r in X n if 
U. - {xtX n ;\\x\\ ^r}. 

The following theorem, by Gohberg and Krein, is proved in Ref. 7 
and will be found useful. 

Theorem: If X n+1 is an (n + l)-dimensional subspace of a Banach 
space (B and if U n+1 is the closed ball of radius r in X„ +1 then d n (U n+1 ) = r. 

This theorem and the nested property of widths can be used to obtain 
lower bounds on d n (5) . This lower bound can be obtained by constructing 
a ball in an (n + l)-dimensional subspace and choosing r such that 
£/ n+ i C ff. Using the nested property then leads to 

r = d„(U n+l ) g «L(ff). (21) 

Since 

d k ^ d k (22) 

the radius of ball also serves as a lower bound on (/*)*. 

Lemma 1: Let $(t, t ) and C(t) be the transition matrix and output 
matrix, respectively, of the original system (1) and (2). Assume this system 
to be completely observable on [t, , t f ). Let W(t) satisfy the previously stated 
conditions. Then the matrix 



M = [" *'(«, t )C'(t)W(t)C(t)$(t, t ) dt 

is positive definite. 

Proof: Consider the quadratic form x' Mx Q = || y || 2 ^ where 
x{t ) = x , that is, y(t) = C{t)$(t, t )x . 
Now II y || 2 = => y(t) ^ on [*. , t,]. 



(23) 
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Since the system is observable y = => x = 0. Thus M is positive 
definite. 

The following theorem provides the lower bound on the performance. 

Theorem 1: Let R be the closed region oj initial conditions on the 
original system and let x(t ) = be an interior point of R. Assume the 
system to be completely observable on [t. , t,\. Denote the boundary of R by 
dR and let 

p 2 4 m in x'(tMQ. (24) 

z(.t°)tdR 

Let the eigenvalues of the positive definite matrix M be ordered Xj(M) ^ 
^(M) ^ • * * s£ X n (M). Then the performance, for any k-dimensional 
approximating system, satisfies J k (A) ^ p 2 \ k+l (M) for k < n. 

Proof: Let 

ff = {y; y(t) = C(t)*(t, 0*(0, *('.) e R). 

A k + 1 dimensional ball will now be constructed which is a subset of 
5\ Consider the k + 1 dimensional ball of radius r 

U k+1 = {y; y(t) = C(t)Ht, 0*(0. *(0 * #* + i C £„, II y II ^ r} 

# i+1 and r will be chosen such that U k+i C &• Since M is real and sym- 
metric it can be diagonalized with an orthogonal matrix T. Thus M = 
T'AT and 

||y|| 2 = [Tx(t.)YA.[Tx(t )] = z'Az 
where 





T' = 






A = 


Lo x n _ 


and 






z = Tx(t ). 


Defining 




Ek+i = 


- \x(t e ); [TxQMt = 0, * 


and 






r~ = 


= P 2 X t+1 (M) 



i = k + 2, 
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gives 

U k+l = [y; y(t) = C(t)$(t, t„)x(Q, 

|| y || 2 ^ pXUM), [7MQ].=0 
i = k -f- 2, — ,n\. 
Thus for y e £/*+, 

|| y || 2 = x'(t )Mx(Q = £s?X< ;g p 2 X 4+1 . 

Since 

z { = * = ft + 2, • • • , n 

and 



X,- 

X*+i 



we have 



^ 1 for i ^ fc + 1 



£z! = z£z ^ p 2 



It then follows, from the definition of p and the fact that zero is an 
interior point of R, that x t R and therefore y e JF. Thus U k+J C & 
and the desired result 

J k (l) }z d 2 k {5) > P \ +1 (M) k<n 

follows. 

Remarks: Recalling that the eigenvalues of M are ordered, we 
notice that the lower bound is a decreasing function of the dimension 
of the approximating system. This result can be used to determine 
what order aproximating system (at least) need be considered to 
achieve a given performance. We emphasize that the bound depends 
on the original system and is obtainable prior to the modelling proce- 
dure. From an engineering viewpoint, if one has an approximating 
system whose performance is "close" to the bound it may not be 
necessary to seek the minor improvement. Notice that the only prop- 
erty of R appearing in the lower bound is p and no attempt was made 
to take the orientation of the set into account. The bound will there- 
fore be least conservative when R is a hypersphere of radius p. 
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IV. EVALUATING THE PERFORMANCE FUNCTION 

In this section the problem of finding the performance (or equi- 
valently the deviation) of a given approximating system is considered. 
The optimal choice of initial conditions on the approximating system 
is obtained using some elementary Hilbert space concepts and it is 
shown that 

inf II y - y II 2 

is a positive semidefinite quadratic form in x(t ). Next, properties of 
convex functions are used to evaluate the performance for different 
classes of regions of initial conditions; namely, for ellipsoids and con- 
vex polyhedra. The Powell algorithm for minimizing a function of 
several variables, without calculating derivatives, is then outlined 
and applied to the system approximation problem. 

The problem of finding 



5 2 = inf || y -y 



(25) 



StVi 



is equivalent to finding the best choice of initial conditions on a given 
approximating system characterized by ^cDi. It can be shown 11 
that there exists a unique y* c of* (y* is called the projection of y in the 
space y k ) such that 

= II y II 2 - II r II 2 - (26) 



tf = 



y - y* 



Furthermore, since g, , g 2 , • • • , g* spans <y* , y* has the representation 

y* = E Bi*f 



where 



G(g. i • ■ • , &)** = 



(y>gO 
_(y, gO. 



(27) 



x* = 
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and G is the Grammian of {g, ;i = 1, • • • , k\, that is, 

[G(g t , ■■■ , g k )]i.t = (g, , gy) i, 3 - li ' " " » k 

Any solution to (27) results in an optimum choice of initial conditions 
on the approximate system. If the g.'s are linearly independent (this 
corresponds to the system being observable) the Grammian is invertible 
and £* is unique. Thus 



**(0 = G\ gl , ■■■ ,g k ) 

where 6 ,f is the pseudoinverse 12 of G. 
The Grammian is given by 



(y, gi) 
L(y, g*)J 



(28) 



Gfei, ■•• ,fc) = f f e I,u - t ' ) C'W(t)Ce Au -' 


' dt 


(29) 


and 






(y, g<) = K'tFziQ 




(30) 


where F is given by 






F = f e*' (t -'' ) C'W(t)C(t)4?(t, to)dt. 




(31) 


Using (30) in (28) gives 






£*{t ) = G*Fx(t ). 




(32) 



Thus the optimal initial condition on the approximating system 
is obtained by linearly transforming the actual initial condition with 
the (k X n) matrix G*F. Using the orthogonality property (26) yields 

l|y-y*ll 2 - \\y\\ 2 -2*'(t )G £ *(t ). 

Letting 

M = f ' $'(£, QC'(t)W(t)C(t)Ht, O dt (33) 

J t. 

and using (32) and the symmetry of G (and thus G f ) gives 

|| y - y* || 2 = x'(Q(M - F'G i F)x(t ). (34) 

In summary, 

5 2 = inf || y - y || 2 = x'{t a ) Dx(Q (35) 

ytH 4 
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with 

D = M - F'G*F. (36) 

Thus, finding the optimal initial condition on the approximating sys- 
tems leads to the positive semidefinite quadratic form (35) for the 
approximation error. The above represents the first step in evaluating 
the performance of any given approximating system. 

Since D is a positive semidefinite matrix, S 2 defined by (35) is a 
convex function of the initial state x(t ). The following theorem from 
Ref.13 is useful in maximizing 8 2 . 

Theorem: If the absolute maximum of a convex function, defined 
on a closed, bounded, convex set, is finite then the absolute maxi- 
mum is taken on at an extreme point of the set. 

Remarks: An extreme point of a convex set is a point in the set 
that cannot be written as a convex combination of two other points 
in the set. Notice that an extreme point is a boundary point; how- 
ever, generally not every boundaiy point is an extreme point. Thus, 
if one is seeking the absolute maximum of a convex function defined 
on a closed, bounded, convex set only boundary points need be con- 
sidered. Also if the domain of definition is a convex polyhedron (a 
closed, bounded, convex set with a finite number of extreme points) 
the absolute maximum can be obtained by simply evaluating the 
function at the extreme points and choosing the largest value. 

Two general classes of closed, bounded, convex regions of initial 
conditions are considered in this paper, the ellipsoid and the convex 
polyhedron. 

Let the region under consideration be an ellipsoid defined by 

R = {x(U; x'(t )Bx{t ) ^ r 2 } (37) 

where B is a positive definite, symmetric matrix and r is finite. Notice 
that R is closed, bounded, and convex. Now the constrained maxi- 
mization problem is one with an inequality constraint. Using the con- 
vexity of R and 8 2 , the absolute maximum of the quadratic form is 
seen to take place on the boundary of the set R. Thus the performance 
can be written 

./,(!) = max x'(Q Dx{t ) 

r(l.) 

subject to the constraint 

x'(t e )Bx(Q = r\ 
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It can easily be shown that the x(t ) maximizing the quadratic form 
is the eigenvector of the matrix B~ X D corresponding to the largest 

eigenvalue and the maximum is given by 

J k (l) = K^iB-'Dy. (38) 

A convex polyhedron is usually representative of the type of in- 
formation one has as to the range of initial conditions. As an example 
of this situation consider the original system to represent linearized 
equations of motion of a space vehicle. Suppose it is known that the 
range of initial conditions are in terms of bounds on position, ve- 
locity deviations, and so on. For example, 

I xM | g 100 feet. 

I x 2 (t ) | ^ 5 feet per second. 

This particular region is described by a rectangular region in state 
space with the extreme points being the corners 



100 




100 




-100 




-100 


5 _ 




_-5. 




. 5 . 




. -5 _ 



In general for this type of initial condition region, that is, 

I Xi(Q | ^ bi i = 1, •• • , n, 

the region has 2" extreme points. Since 8 2 is an even function of x(t ) 
it is only necessary to consider 2"- 1 extreme points eliminating from 
consideration the negative of any point considered. 

The convex polyhedron region also is important, for example, since 
it may be used to simply approximate a more complex region. In gen- 
eral, let 



.(«') 



i = 1, 2, 



,N 



be the extreme points of the convex polyhedron R. Using the con- 
vexity of 8 2 in the initial state x(t ) the absolute maximum 8 2 over 
R takes place at one of the x {i) . Letting 

5* = x (ii 'Dx li) 
where D is given by equation (36) leads to 

J k (l) = maxLtf , hi , ••• , A], (39) 
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V. MINIMIZING THE PERFORMANCE FUNCTION 

Since it is a fairly simple matter to evaluate the performance, whereas 
evaluating the gradient of the performance function requires significant 
computational effort, it is desirable to use a numerical procedure not 
requiring a gradient computation. Notice that Jk(A) is not generally 
differentiable. Here, for completeness, the Powell method of minimizing 
a function of several variables without calculating derivatives is pre- 
sented. 14 Reference 15 contains a summary of the various minimization 
techniques available not requiring the computation of a derivative. 
See Refs. 14 and 15 for a more detailed description of the methods 
and their convergence properties. 

Consider a real, scalar, valued function of N real variables a^ , • ■ ■ , a N 
written /(a). Powell's iterative scheme concerns itself with finding the 
minimum of /(a) without computing its derivative. 

Each iteration of the modified Powell procedure starts with a search 
down N linearly independent directions 

starting with an initial guess a and defines a new set of directions for 
the next iteration. 
An iteration of the recommended procedure, suggested by Powell, is: 

(i) for j = 1, 2, ■ ■ • , N calculate X; such that /(a,_i -f- \flt) is 
minimum and define a t = aj_i + X;ij/ . 

(w) Find the integer m, 1 ^ ra ^ N, such that /(a m _,) — /(a m ) is 
a maximum and define A = /(a m _j) — /(a m ). 
(Hi) Calculate f 3 = f(2a N — a ) and define 

U - /(a.) 
U - 1(flir)> 
(iv) If either / 3 ^ j x or 

(/x - 2/ 2 + /,)(/, - U ~ A) 2 ^ *A(/, - U? 

use the old directions i)x , • • • , t] N for the next iteration and use a N 
for the next a„ , otherwise 

(v) define t\ = a„ — a„ and calculate X such that j(a N + Xrj) is 
minimum. Use 

as the new directions and a N + X17 as the new a . 
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The performance functions, for the two classes of initial conditions 
being considered are given by (38) and (39) in terms of the matrix D 
defined in (36). The major effort in computing the performance func- 
tion is seen to he in the computation of D. Sylvester's expansion (see 
page 83 of Ref. 16) for computing e Xt is useful in the computation of 
the matrices F and G. 

The basic procedure can be outlined as follows: 

(i) Compute and store C(t)$(t, t„) for t t [t, , t f ] using (7) and (8). 
(ii) Evaluate M using (23). 

(Hi) If it is desired to compute the lower bounds to aid in choosing 
the dimension of the approximating system, compute the eigenvalues 
of M and obtain the bounds from the result of Theorem 1. 

(iv) Choose starting values for A and choose the directions for the 
initial search in the modified Powell method to be 



[~l] 




r °l 




r°i 





1 


1 




» 1 





_0_ 




_0_ 




,1_ 



where the above are fc 2 vectors. 

(v) Use modified Powell method to determine the minimum of the 
performance function. Each element of the vector a in the Powell 
method corresponds to an element of A. 

VI. EXAMPLE 

A linearized missile guidance loop may be expressed in the form 



Xl U/2 y 



x 2 = 



H 



m — t 



x z , 



x 3 = u 



(40) 



where Xi is the lateral position deviation from a nominal trajectory, 
x 2 is the lateral velocity deviation, x 3 is the attitude deviation in the 
given direction and u is the control signal. The relationship between 
the attitude and lateral acceleration is given through the time-varying 
gain H/(m — t) which accounts for the loss of mass because of fuel 
consumption. 



MODELING CONTROL SYSTEMS 



225 



Suppose it is desired to approximate homogeneous solutions to 
(40) for initial conditions (at beginning of a stage) lying in a set 
R (R is denned later) with solutions of a constant coefficient system. 
The actual system (40) can be written in the vector-matrix form 



with output 



where 



x(t) = A(t)x(t) 
y(t) = [10 0}x(t) = Cx(t) 



(41) 
(42) 



x{t) = 



x 2 (t) 

Le,(l)J 



and 



A(t) = 



r° 


1 











H 

m — t 


_0 





_ 



(43) 



Let 



y || 2 = [ T y 2 (t)dt. 

Jo 



(44) 



Before proceeding to find the approximation it is instructive to 
determine the lower bounds on the optimal performance. This will 
naturally aid in choosing the dimension of the approximating system. 
The matrix, M, defined by (33) , is given by 



M 



= ( t>'(t, o)C'C${t, o) 



dt 



with 



d_ 

dt 



*(t, o) = A(t)*(t, o). 



(45) 



(46) 



The transition matrix, which is the solution to (46) with the identity 
initial condition, is given by 
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Evaluating M leads to 



1 t HUm - Olnl^r^j -!- / 



1 
.0 



-tf In 



m 

m — t 
m 



M = 



T £ *£-^l» ^) + i ( r- 2m T) 



/T72 /Tj3 

y y 

Lm 13 m 2 



M 23 
M 33 



with 



Ma ^ g r£_5 3 _ (2r + m)(m- / 



36 



•In 



6 



w - T\ . (m - T) 2 (4T + 5m) 



to 



+ 



36 



and 



M„ = tf 



!F 3 (m - T) 3 ., (m-T 



m 



lnM ^ ±.j + ^ (m _ jp). 



(2T + to) 



10 3 
~36 m " 



(to - D s 



! ! m ~ T ) . ( m ~ :r ) 2 ( 471 + 5w ) 



??i 



18 



Let the constants defining the problem be given by 

wi = 15 seconds (normalized mass) 

T = 10 seconds 

i/ = 15 (pound-seconds per slug ) X 10 -3 
and let the region of initial conditions be given by 

R = \x{o); | x x (o) | ^ 30 feet, | x 2 (o) | ^ 2 feet per second, 

I £3(0) I ^ 1 milliradian } . 
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Evaluating M for the above values of the constants leads to 

10 50 206" 

1/ 50 333 1570 

206 1570 8082 
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with eigenvalues 



We have 



and 



X, = 8393, X 2 = 31, X 3 = 1.1. 

J ^ 8,393 
Jx ^ 31 

J, ^ 1.1. 



Here Jo represents 



max 1 1 y 



The second order approximation thus has the possibility of yielding 
a negligible approximation error. Thus in the remainder of this paper 
the optimal second order approximation will be sought. Thus 



£ = Ax = 



(In Q, l2 

l_t*21 **2 2 _] 



and 



y = [l 0}£. 

The initial choice for A in the iterative procedure is 

l" 



I. - 







which represents polynomial approximations to solutions of the orig- 
inal system. 
The extreme points of R are given by 





30 




-30 




30 






30 


x (,) - 


2 
1 


, * W = 


•2 
1_ 


, * <3) = 


-2 
1 


and 


x w = 


2 
-1 
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and their negatives. Thus 

J a (I„) = max [x w ' Dx w ) = 340 
i 

where D is evaluated from (36) . It is thus seen that the performance 
function is far greater than the lower bound and the possibility exists 
for a significant improvement. The result of applying the Powell 
algorithm to this problem yields 



A* = 



0.244 0.827 

0.177 X 10" 8 0.629 X 10" 8 J 



and 



J 2 (A*) = 33.4 
with the eigenvalues of A* given by 

X x (l*) = 0.245 

X 2 (2*) = 0.30 X 10~ 4 . 

The above results are obtained after three iterations of the Powell 
algorithm. The G, F and D matrices are given by 

"271 77l" 

.771 2230. 

"43.14 296.0 1459" 

_112.2 832.6 4241. 

and 



G = 



F = 





4.4 X 10~ fl 


-1.8 X 10~* 


1.2 X 10" 4 


D = 


-1.8 X 10" 4 


7.3 


3.5 X HT 3 




1.2 X 10" 4 


3.5 X 10" 1 


4.2 


Evaluating 










max 
i 


[x tn Dx li) ) 





gives the maximum approximation error occurring at the extreme point 



30 

-2 

1_ 
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Figure 1 shows the solution of the actual and approximate system 
for this worst-case initial condition. The solutions are obtained from 



y(t) = 30 - 2t + *i 3 (i, o) 



and 



I««/-!-li 



m = [1 0V'G- l F 



30 

-2 

1 



y(t) = 4.82 e x '' + 19.78 e x,t . 
The matrix relating the initial conditions is given by G~ r F, that is, 

£(o) = G~ l Fx{o). 

1.00 1.86 -1.68" 1 

-0.295 -0.271 2. 48 J 



£(o) = 



x(o). 




Fig. 1 — Exact and approximate solutions in worst case. 



230 THE BELL SYSTEM TECHNICAL JOURNAL, JANUARY 1969 

vn. CONCLUSIONS 

A method for uniformly approximating solutions of linear, time- 
varying, homogeneous differential equations has been presented. The 
problem of approximating systems subject to control or reference 
inputs is considered in Ref. 6 for the class of exponential polynomial 
control inputs. 

One of the objectives of modeling with constant coefficient systems 
was to obtain closed form approximations. Use of Sylvester's expan- 
sion allows one to derive these closed form expressions. However, 
more general classes of approximating systems can be sought while 
still maintaining the property that approximations are in closed form. 
For example, a general model of the form 

£ = p(t)A£ 

where p{t) is a scalar valued function possesses the closed form 
solution 



[l f p( T ) dr]s 



£{t) = exp \A / p(r) dr mt ) 

and p{l) as well as A may be sought as part of the modeling procedure. 
A complexity constraint can be imposed on p(t) by considering it to 
be a polynomial of given degree and the search for the model reduces 
to finding the coefficients of the polynomial as well as the elements of A . 
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